import numpy as np
from scipy.stats import norm

# ================= 参数设置 =================
S0 = 100      # 初始股价
K = 100       # 执行价
r = 0.05      # 无风险利率
sigma = 0.2   # 波动率
T = 1         # 到期时间（年）
N = 1000      # 模拟路径数
# np.random.seed(42)  # 保证结果可复现

# ================= 1. 蒙特卡洛模拟估计期权价格 C_MC =================
# 生成标准正态随机数
Z = np.random.standard_normal(N)

# 几何布朗运动模拟到期股价 S_T
ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)

# 计算到期收益（看涨期权）
VT = np.maximum(ST - K, 0)

# 蒙特卡洛估计的期权价格
C_MC = np.exp(-r * T) * np.mean(VT)
print(f'蒙特卡洛估计价格 C_MC: {C_MC:.4f}')

# ================= 2. Black-Scholes 理论价格 C_BS 和相对误差 =================
def black_scholes_call(S0, K, r, sigma, T):
    d1 = (np.log(S0/K) + (r + 0.5*sigma**2)*T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    C = S0 * norm.cdf(d1) - K * np.exp(-r*T) * norm.cdf(d2)
    return C, d1, d2

C_BS, d1, d2 = black_scholes_call(S0, K, r, sigma, T)
relative_error = abs(C_MC - C_BS) / C_BS

print(f'BS 理论价格 C_BS: {C_BS:.4f}')
print(f'相对误差: {relative_error:.4f}')

# ================= 3. 控制变量法 (Control Variates) =================
# 控制变量：到期股价 S_T，其理论期望为 S0 * exp(r*T)
E_ST = S0 * np.exp(r * T)  # 理论期望
sim_ST_mean = np.mean(ST)  # 模拟均值

# 计算协方差 Cov(ST, VT) 和方差 Var(ST)
cov_ST_VT = np.cov(ST, VT)[0, 1]  # 协方差
var_ST = np.var(ST)            # ST 的方差

# 最优控制变量系数 beta
alpha = cov_ST_VT / var_ST

# 控制变量法修正后的期权价格
C_CV = C_MC - alpha * (sim_ST_mean - E_ST)

print(f'\n控制变量法估计价格 C_CV: {C_CV:.4f}')
print(f'控制变量法相对误差: {abs(C_CV - C_BS) / C_BS:.4f}')

# 计算两种方法的样本方差（评估效率）
var_MC = np.var(np.exp(-r*T) * VT) / N
var_CV = np.var(np.exp(-r*T) * (VT - alpha * (ST - E_ST))) / N

print(f'\n方差比较（每条路径）:')
print(f'  标准 MC 方差: {var_MC * N:.4f}')
print(f'  控制变量法方差: {var_CV * N:.4f}')
print(f'方差减少比例: {1 - var_CV / var_MC:.2%}')

